module Index

use strutils
use OptimisationTypes
use SysUtils 


implicit none


contains

function GetIndexFileName(Opt)
   implicit none
   type(TOptimisation) :: Opt
   character(len = 256) :: GetIndexFileName
   GetIndexFileName = GetWorkFileName(Opt, s_IndexFile)
end function GetIndexFileName


function GetDataKeyFileName(Opt)
   implicit none
   type(TOptimisation) :: Opt
   character(len = 256) :: GetDataKeyFileName
   GetDataKeyFileName = GetDataFileName(Opt, s_DataKeyFile)
end function GetDataKeyFileName


subroutine CreateNewIndexData(Opt, FirstYear, LastYear, FirstAge, LastAge, NumFleets, NumSSB, NumBiomass)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iDataSourceType, NumElms, I
   integer :: FirstYear, LastYear, FirstAge, LastAge, NumFleets, NumSSB, NumBiomass
   type(TDataSourceArray), pointer :: pDataSourceArray
   type(TDataSource), pointer :: pDataSource
   
   ! Reallocate index data
   call DeallocateIndexData(Opt)
   allocate(Opt%IndexInput)
   allocate(Opt%IndexInput%DataSources(0 : size(c_DataSource) - 1))

	Opt%IndexInput%DataFiles(odCaton) = "caton.txt"
	Opt%IndexInput%DataFiles(odCanum) = "canum.txt"
	Opt%IndexInput%DataFiles(odWeca) = "weca.txt"
	Opt%IndexInput%DataFiles(odWest) = "west.txt"
	Opt%IndexInput%DataFiles(odMatProp) = "matprop.txt"
	Opt%IndexInput%DataFiles(odFleet) = "fleet.txt"
	Opt%IndexInput%DataFiles(odSSB) = "ssb.txt"
	Opt%IndexInput%DataFiles(odBiomass) = "biomass.txt"
   Opt%IndexInput%StockTitle = "Stock assessment"
   Opt%IndexInput%ProcessType = ptOrdinaryAssessment
   Opt%IndexInput%SearchRoutine = srGeneric
   Opt%IndexInput%ModelType = pmtVPA
   Opt%IndexInput%ISVPAControlType = ISVPA_CatchControl
   Opt%IndexInput%PlusGroupMethod = pgDynamicPool
   Opt%IndexInput%FirstYear = FirstYear
   Opt%IndexInput%LastYear = LastYear
   Opt%IndexInput%LastCatchYear = LastYear
   Opt%IndexInput%FirstAge = FirstAge
   Opt%IndexInput%LastAge = LastAge
   Opt%IndexInput%FirstRefAge = FirstAge
   Opt%IndexInput%LastRefAge = LastAGe
   Opt%IndexInput%NumMedianElms = 10
   !Opt%IndexInput%GammaShape = 0.0d0
   Opt%IndexInput%NumStages = 1000
   Opt%IndexInput%SSBFracF = 0.0d0;
   Opt%IndexInput%SSBFracM = 0.0d0;
   do iDataSourceType = 0, size(c_DataSource) - 1
      pDataSourceArray => Opt%IndexInput%DataSources(iDataSourceType)
      select case(c_DataSource(iDataSourceType))
         case(odCanum) 
            NumElms = 1
         case(odFleet) 
            NumElms = NumFleets
         case(odSSB)            
            NumElms = NumSSB
         case(odBiomass)
            NumElms = NumBiomass
      end select
      if (NumElms == 0) cycle
      allocate(pDataSourceArray%Elm(0 : NumElms - 1)) 
      do I = 0, NumElms - 1
         pDataSource => pDataSourceArray%Elm(I)
         pDataSource%FProp = 0.0d0
         pDataSource%MProp = 0.0d0
         pDataSource%ObjType = otSSQ
         pDataSource%Weight = 1.0d0
         pDataSource%My = 0.001d0
         pDataSource%PlusGroupCode = 1
            pDataSource%BootType = btOrgRes
         pDataSource%Sigma = 0.0d0
         pDataSource%BootFile = "NA"
      enddo
   enddo
   Opt%IndexInput%InParamFileName = "paramin.txt"
   Opt%IndexInput%DataFiles(odCanumWeight) = "CanumWgt.txt"
   Opt%IndexInput%DataFiles(odFleetWeight) = "FleetWgt.txt"
   Opt%IndexInput%DataFiles(odSSBWeight) = "SSBWgt.txt"
   Opt%IndexInput%DataFiles(odBiomassWeight) = "BiomWgt.txt"

   Opt%IndexInput%OutFiles(roParamOut) = "paramOut.txt"       
   Opt%IndexInput%OutFiles(roSummary) = "Summary.txt"       
   Opt%IndexInput%OutFiles(roOutput) = "Output.txt"       
   Opt%IndexInput%OutFiles(roBootResult) = "Boot_Result.txt"       
   Opt%IndexInput%OutFiles(roBootFMean) = "Boot_FMean.txt"       
   Opt%IndexInput%OutFiles(roBootSSB) = "Boot_SSB.txt"       
   Opt%IndexInput%OutFiles(roBootObj) = "Boot_ObjFunc.txt"       
   Opt%IndexInput%OutFiles(roRetroFMean) = "Retro_FMean.txt"       
   Opt%IndexInput%OutFiles(roRetroSSB) = "Retro_SSB.txt"       
   Opt%IndexInput%OutFiles(roRetroRecruit) = "Retro_Recruit.txt"       
end subroutine CreateNewIndexData


function ReadIndexData(Opt)
! Purpose: Read input data from index file
!          Index file name is entered by user
! History: When       Who       What
!          13.04.07   aasmund   Created
   implicit none
   logical :: ReadIndexData
   type(TOptimisation), intent(inout) :: Opt
   integer :: ioerr, i, iDataSourceType, iDataSource, NumElms, s
   type(TDataSourceArray), pointer :: pDataSourceArray
   type(TDataSource), pointer :: pDataSource
   character(len = 80) :: FirstLine
   double precision :: GammaShape
   ReadIndexData = .false.
   ioerr = -1

   ! Reallocate index data
   call DeallocateIndexData(Opt)
   allocate(Opt%IndexInput)
   allocate(Opt%IndexInput%DataSources(0 : size(c_DataSource) - 1))

   if(.not. OpenFile(Opt, cIndexFile, GetIndexFileName(Opt), .false.)) return
   !open(cIndexFile, file = GetIndexFileName(Opt), action = 'read', iostat = ioerr, status = 'old')
   ! Read the stock title
	
	! Read the main line
   ! Read Model variables
   if(.not. ReadStrFromFile(Opt, cIndexFile, FirstLine)) return
   if(FirstLine(1:3) == 'VER') then
      Opt%FileVersion = strtofloat(FirstLine(5:7))
      if(.not. ReadStrFromFile(Opt, cIndexFile, Opt%IndexInput%StockTitle)) return
   else
      Opt%FileVersion = 1.0
      Opt%IndexInput%StockTitle = FirstLine
   endif

	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%ProcessType, 0, 3, ElmName = "Process type")) return
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%SearchRoutine, 0, 0, "Search routine")) return
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%ModelType, 0, 2, "Model type")) return
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%ISVPAControlType, 0, 1, "ISVPA type")) return
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%PlusGroupMethod, 0, 1, "Plus group method", .true.)) return

	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%FirstYear, 0, -1, "First year")) return
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%LastYear, 0, -1, "Last year")) return
   if(.not. CheckRange(Opt, Opt%IndexInput%FirstYear, Opt%IndexInput%LastYear, "First year", "Last year")) return
	
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%LastCatchYear, &
	   Opt%IndexInput%FirstYear, Opt%IndexInput%LastYear, "Last catch year")) return
	
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%FirstAge, 0, -1, "First age")) return
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%LastAge, 0, -1, "Last age")) return
   if(.not. CheckRange(Opt, Opt%IndexInput%FirstAge, Opt%IndexInput%LastAge, "First age", "Last age")) return
	
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%FirstRefAge, &
	   Opt%IndexInput%FirstAge, Opt%IndexInput%LastAge, "First reference age")) return
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%LastRefAge, &
	   Opt%IndexInput%FirstAge, Opt%IndexInput%LastAge, "Last reference age", .true.)) return
   if(.not. CheckRange(Opt, Opt%IndexInput%FirstRefAge, Opt%IndexInput%LastRefAge, "First reference age", &
   	"Last reference age")) return

   if(Opt%FileVersion >= 1.3) then
      if(.not. ReadFFromFile(Opt, cIndexFile, Opt%IndexInput%SSBFracF, 0.0d0, -1.0d0, "SSBFracF")) return
      if(.not. ReadFFromFile(Opt, cIndexFile, Opt%IndexInput%SSBFracM, 0.0d0, -1.0d0, "SSBFracM", .true.)) return
   else
      Opt%IndexInput%SSBFracF = 0.0d0;
      Opt%IndexInput%SSBFracM = 0.0d0;
   endif
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%NumMedianElms, 0, -1, "Number of median elements")) return
	if(Opt%FileVersion < 1.2d0) then
	   if(.not. ReadFFromFile(Opt, cIndexFile, GammaShape, -1.0d0, -1.0d0, "Gamma shape")) return
	endif
	if(.not. ReadIFromFile(Opt, cIndexFile, Opt%IndexInput%NumStages, 1, -1, "Number of process stages", .true.)) return
   do iDataSourceType = 0, size(c_DataSource) - 1
      pDataSourceArray => Opt%IndexInput%DataSources(iDataSourceType)
   	if(.not. ReadIFromFile(Opt, cIndexFile, NumElms, 0, 100, "Number of datasources", .true.)) return
      if (NumElms == 0) cycle 
      allocate(pDataSourceArray%Elm(0 : NumElms - 1)) 
      do iDataSource = 0, size(pDataSourceArray%Elm) - 1  
         pDataSource => pDataSourceArray%Elm(iDataSource)
      	if(.not. ReadFFromFile(Opt, cIndexFile, pDataSource%FProp, 0.0d0, 1.0d0, "Proportion of F")) return
      	if(.not. ReadFFromFile(Opt, cIndexFile, pDataSource%MProp, 0.0d0, 1.0d0, "Proportion of M")) return
      	if(.not. ReadIFromFile(Opt, cIndexFile, pDataSource%ObjType, 0, 5, "Objective function type")) return
      	if(.not. ReadFFromFile(Opt, cIndexFile, pDataSource%Weight, 0.0d0, 1000d0, "Datasource weight")) return
      	if(Opt%FileVersion == 1.0) then
      	   if((iDataSourceType == 0) .and. (Opt%IndexInput%ModelType /= pmtSep)) then
      	      ! The old VPA/ISVPA catch residuals was allways ignored in Tasacs 1.0
      	      ! 1.1 Compatibility to this is to use datasource weight = 0
      	      pDataSource%Weight = 0.0d0
      	   endif
      	endif
      	if(.not. ReadFFromFile(Opt, cIndexFile, pDataSource%My, 0.0d0, 1.0d0, "Objective my constant")) return
      	if(.not. ReadIFromFile(Opt, cIndexFile, pDataSource%PlusGroupCode, 0, 1, "Plus group code" )) return
      	if(.not. ReadIFromFile(Opt, cIndexFile, pDataSource%BootType, 0, 4, "Bootstrap type")) return
      	if(.not. ReadFFromFile(Opt, cIndexFile, pDataSource%Sigma, 0.0d0, 1.0d0, "Sigma")) return
      	if(.not. ReadStrFromFile(Opt, cIndexFile, pDataSource%BootFile)) return
      enddo
   enddo
   ! Run input files
   if(.not. ReadStrFromFile(Opt, cIndexFile, Opt%IndexInput%InParamFileName)) return
   if(.not. ReadStrFromFile(Opt, cIndexFile, Opt%IndexInput%DataFiles(odCanumWeight))) return
   if(.not. ReadStrFromFile(Opt, cIndexFile, Opt%IndexInput%DataFiles(odFleetWeight))) return
   if(.not. ReadStrFromFile(Opt, cIndexFile, Opt%IndexInput%DataFiles(odSSBWeight))) return
   if(.not. ReadStrFromFile(Opt, cIndexFile, Opt%IndexInput%DataFiles(odBiomassWeight))) return
   ! Run output files
   do i = roFirst, roLast
      if(.not. ReadStrFromFile(Opt, cIndexFile, Opt%IndexInput%OutFiles(i))) return
	end do
   close(cIndexFile)
   
   ! Read data key file
   if(.not. OpenFile(Opt, cDataKeyFile, GetDataKeyFileName(Opt), .false.)) return
	if(.not. ReadStrFromFile(Opt, cDataKeyFile, Opt%IndexInput%DataFiles(odCaton))) return
	if(.not. ReadStrFromFile(Opt, cDataKeyFile, Opt%IndexInput%DataFiles(odCanum))) return
	if(.not. ReadStrFromFile(Opt, cDataKeyFile, Opt%IndexInput%DataFiles(odWeca))) return
	if(.not. ReadStrFromFile(Opt, cDataKeyFile, Opt%IndexInput%DataFiles(odWest))) return
	if(.not. ReadStrFromFile(Opt, cDataKeyFile, Opt%IndexInput%DataFiles(odMatProp))) return
	if(.not. ReadStrFromFile(Opt, cDataKeyFile, Opt%IndexInput%DataFiles(odFleet))) return
	if(.not. ReadStrFromFile(Opt, cDataKeyFile, Opt%IndexInput%DataFiles(odSSB))) return
	if(.not. ReadStrFromFile(Opt, cDataKeyFile, Opt%IndexInput%DataFiles(odBiomass))) return
   close(cDataKeyFile)

   ReadIndexData = .true.
end function ReadIndexData


subroutine WriteIndexData(Opt)
! Purpose: Write input data to index file
! History: When       Who       What
!          30.05.07   aasmund   Created
   implicit none
   type(TOptimisation), intent(in) :: Opt
   integer :: ioerr, i, iDataSourceType, iDataSource, NumElms
   type(TDataSourceArray), pointer :: pDataSourceArray
   type(TDataSource), pointer :: pDataSource

   ioerr = -1
   open(cIndexFile, file = GetIndexFileName(Opt), action = 'write', iostat = ioerr, status = 'replace')
   if (io_no_error .eq. ioerr) then
      ! Write the version number
      write(cIndexFile, '(A)') s_CurrentVersion
      ! Write the stock title
      write(cIndexFile, '(A)') Opt%IndexInput%StockTitle
      ! Write Model variables
      write(cIndexFile, "(5I10)") &
         Opt%IndexInput%ProcessType, &
         Opt%IndexInput%SearchRoutine, &
         Opt%IndexInput%ModelType, &
         Opt%IndexInput%ISVPAControlType, &
         Opt%IndexInput%PlusGroupMethod
      write(cIndexFile, "(7I10)") &
         Opt%IndexInput%FirstYear, &
         Opt%IndexInput%LastYear, &
         Opt%IndexInput%LastCatchYear, &
         Opt%IndexInput%FirstAge, & 
         Opt%IndexInput%LastAge, &
         Opt%IndexInput%FirstRefAge, &
         Opt%IndexInput%LastRefAge
      write(cIndexFile, "(2F10.3)") & 
         Opt%IndexInput%SSBFracF, &
         Opt%IndexInput%SSBFracM
      write(cIndexFile, "(I10I10)") &
         Opt%IndexInput%NumMedianElms, & 
         !Opt%IndexInput%GammaShape, &
         Opt%IndexInput%NumStages
      do iDataSourceType = 0, size(c_DataSource) - 1
         pDataSourceArray => Opt%IndexInput%DataSources(iDataSourceType)
         NumElms = 0
         if(associated(pDataSourceArray%Elm)) then
            NumElms = size(pDataSourceArray%Elm)
         endif
         write(cIndexFile, "(I10)") NumElms
         if (NumElms == 0) cycle
         do iDataSource = 0, size(pDataSourceArray%Elm) - 1  
            pDataSource => pDataSourceArray%Elm(iDataSource)
            write(cIndexFile, "(2F10.3I5F10.3F10.5I5I5F10.3A32)") & !Opt%IndexInput%FleetIndices(i)%Weight, &
               pDataSource%FProp, pDataSource%MProp, &
               pDataSource%ObjType, pDataSource%Weight, &
               pDataSource%My, pDataSource%PlusGroupCode, &
               pDataSource%BootType, pDataSource%Sigma, adjustr(pDataSource%BootFile)
         enddo
      enddo
      ! Run input files
      write(cIndexFile, '(A)') Opt%IndexInput%InParamFileName
      write(cIndexFile, '(A)') Opt%IndexInput%DataFiles(odCanumWeight)
      write(cIndexFile, '(A)') Opt%IndexInput%DataFiles(odFleetWeight)
      write(cIndexFile, '(A)') Opt%IndexInput%DataFiles(odSSBWeight)
      write(cIndexFile, '(A)') Opt%IndexInput%DataFiles(odBiomassWeight)
      ! Run output files
      do i = roFirst, roLast
			write(cIndexFile, '(A)') Opt%IndexInput%OutFiles(i)
		end do
      close(cIndexFile)
   else
      call MessageOut(Opt, "Error reading: " // GetIndexFileName(Opt))
   endif

   ioerr = -1
   open(cDataKeyFile, file = GetDataKeyFileName(Opt), action = 'write', iostat = ioerr, status = 'replace')
	write(cDataKeyFile, '(A)') Opt%IndexInput%DataFiles(odCaton)
	write(cDataKeyFile, '(A)') Opt%IndexInput%DataFiles(odCanum)
	write(cDataKeyFile, '(A)') Opt%IndexInput%DataFiles(odWeca)
	write(cDataKeyFile, '(A)') Opt%IndexInput%DataFiles(odWest)
	write(cDataKeyFile, '(A)') Opt%IndexInput%DataFiles(odMatProp)
	write(cDataKeyFile, '(A)') Opt%IndexInput%DataFiles(odFleet)
	write(cDataKeyFile, '(A)') Opt%IndexInput%DataFiles(odSSB)
	write(cDataKeyFile, '(A)') Opt%IndexInput%DataFiles(odBiomass)
   close(cDataKeyFile)

end subroutine WriteIndexData


function FirstIndex() 
! Purpose: Returns the first index.
	implicit none
   integer :: FirstIndex
   FirstIndex = 1
end function FirstIndex


function NumYears(IndexInput) 
! Purpose: Returns the number of years
	implicit none
   type(TIndexInput), intent(in) :: IndexInput
   integer :: NumYears
   NumYears = IndexInput%LastYear - IndexInput%FirstYear + 1;
end function NumYears


function NumCatchYears(IndexInput) 
! Purpose: Returns the number of catch years
	implicit none
   type(TIndexInput), intent(in) :: IndexInput
   integer :: NumCatchYears
   NumCatchYears = IndexInput%LastCatchYear - IndexInput%FirstYear + 1;
end function NumCatchYears


function NumBackwardProjectedYears(IndexInput) 
! Purpose: Returns the number of years to be backward calculated
	implicit none
   type(TIndexInput), intent(in) :: IndexInput
   integer :: NumBackwardProjectedYears
   select case(IndexInput%ModelType)
      case(pmtVPA, pmtISVPA)
         NumBackwardProjectedYears = NumCatchYears(IndexInput)
      case default
         NumBackwardProjectedYears = NumYears(IndexInput)
   end select      
end function NumBackwardProjectedYears


function NumTrueAges(IndexInput)
! Purpose: Returns the number of true ages
	implicit none
   type(TIndexInput), intent(in) :: IndexInput
   integer :: NumTrueAges
   NumTrueAges = NumAges(IndexInput) - 1 ! Minus plus group
end function NumTrueAges


function NumAges(IndexInput)
! Purpose: Returns the number of ages (includes the age plus group)
	implicit none
   type(TIndexInput), intent(in) :: IndexInput
   integer :: NumAges
   NumAges = IndexInput%LastAge - IndexInput%FirstAge + 1
end function NumAges


function NumFleetIndices(IndexInput) 
! Purpose: Returns the number of fleet indices 
	implicit none
   type(TIndexInput), intent(in) :: IndexInput
   integer :: NumFleetIndices
   NumFleetIndices = 0
   if(associated(IndexInput%DataSources(1)%Elm)) then
      NumFleetIndices = size(IndexInput%DataSources(1)%Elm)
   endif
end function NumFleetIndices 


function NumSSBIndices(IndexInput)
! Purpose: Returns the number of SSB indices 
	implicit none
   type(TIndexInput), intent(in) :: IndexInput
   integer :: NumSSBIndices
   NumSSBIndices = 0
   if(associated(IndexInput%DataSources(2)%Elm)) then
   	NumSSBIndices = size(IndexInput%DataSources(2)%Elm)
   endif
end function NumSSBIndices 


function NumBiomassIndices(IndexInput) 
! Purpose: Returns the number of biomass indices 
	implicit none
   type(TIndexInput), intent(in) :: IndexInput
   integer :: NumBiomassIndices
   NumBiomassIndices = 0
   if(associated(IndexInput%DataSources(3)%Elm)) then
      NumBiomassIndices = size(IndexInput%DataSources(3)%Elm)
   endif
end function NumBiomassIndices


function GetDataSource(Opt, DataSourceType, DataSourceIndex)
   implicit none
   type(TOptimisation), intent(in) :: Opt
   integer, intent(in) :: DataSourceType, DataSourceIndex
   integer :: iDataSourceType
   type(TDataSource), pointer :: GetDataSource
   do iDataSourceType = 0, size(c_DataSource) - 1 
      if (c_DataSource(iDataSourceType) == DataSourceType) then
         if (.not. associated(Opt%IndexInput%DataSources(iDataSourceType)%Elm)) exit
         if (size(Opt%IndexInput%DataSources(iDataSourceType)%Elm) == 0) exit
         GetDataSource => Opt%IndexInput%DataSources(iDataSourceType)%Elm(DataSourceIndex)
         return
      endif
   enddo
   GetDataSource => null()
end function GetDataSource


subroutine DeallocateIndexData(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iDataSourceType, NumDataSources, iDataSource
	type(TDataSource), pointer :: pDataSource
	type(TDataSourceArray), pointer :: pDataSourceArray

   if (associated(Opt%IndexInput)) then
      do iDataSourceType = 0, size(c_DataSource) - 1
         pDataSourceArray => Opt%IndexInput%DataSources(iDataSourceType)
         if (associated(pDataSourceArray)) then
            if (associated(pDataSourceArray%Elm)) then
               NumDataSources = size(pDataSourceArray%Elm)
               do iDataSource = 0, NumDataSources - 1
                  pDataSource => pDataSourceArray%Elm(iDataSource)
                  ! Clean up
                  if (associated(pDataSource%BootObs)) then
                     deallocate(pDataSource%BootObs) 
                  endif
               enddo
               deallocate(pDataSourceArray%Elm)
            endif
         endif
      enddo
      if (associated(Opt%IndexInput%DataSources)) then
         deallocate(Opt%IndexInput%DataSources)
      endif
      deallocate(Opt%IndexInput)
   endif
end subroutine DeallocateIndexData


end module Index
